#include <xmmintrin.h> //SSE
#include <emmintrin.h> //SSE2
#include <pmmintrin.h> //SSE3
#include <tmmintrin.h> //SSSE3
#include <smmintrin.h> //SSE4.1
#include <nmmintrin.h> //SSSE4.2
#include <immintrin.h> //AVX?1?7?1?7AVX2
#include <pthread.h>
#include <omp.h>
#include <iostream>
#include <cmath>
#include <semaphore.h>
#include <stdio.h>
#include <windows.h>
#define ROW 1024
#define TASK 10
#define INTERVAL 10
using namespace std;
float matrix[ROW][ROW];
typedef long long ll;
typedef struct {
	int k;
	int t_id;
}threadParam_t;

sem_t sem_main;
sem_t sem_workerstart;
sem_t sem_workerend;
sem_t sem_leader;
sem_t sem_Divsion[32];
sem_t sem_Elimination[32];
pthread_barrier_t division;
pthread_barrier_t elemation;
int NUM_THREADS = 8;
int remain = ROW;
pthread_mutex_t remainLock;

void init()
{
	for(int i = 0;i < ROW;i++)
	{
		for(int j = 0;j < i;j++)
        {
            matrix[i][j] = 0;
        }
		for(int j = i;j<ROW;j++)
        {
            matrix[i][j] = rand() % 10000 + 1;
        }
	}
	for(int i = 0;i < 2000;i++)
	{
		int row1 = rand() % ROW;
		int row2 = rand() % ROW;
		int judge = rand() % 2;
		if(judge == 1)
        {
            for(int j = 0;j < ROW;j++)
            {
                matrix[row1][j] += matrix[row2][j] * (rand() % 100);
            }
        }
		else
        {
            for(int j = 0;j < ROW;j++)
            {
                matrix[row1][j] -= matrix[row2][j] * (rand() % 100);
            }
        }
	}
}

void plain() {
	for (int i = 0; i < ROW - 1; i++) {
		for (int j = i + 1; j < ROW; j++) {
			matrix[i][j] = matrix[i][j] / matrix[i][i];
		}
		matrix[i][i] = 1;
		for (int k = i + 1; k < ROW; k++) {
			for (int j = i + 1; j < ROW; j++) {
				matrix[k][j] = matrix[k][j] - matrix[i][j] * matrix[k][i];
			}
			matrix[k][i] = 0;
		}
	}
}

void SIMD()
{
    for(int k = 0; k < ROW; ++k)
	{
		__m128 diver = _mm_load_ps1(&matrix[k][k]);
		int j;
		for (j = k + 1;j < ROW && ((ROW - j) & 3);++j)
        {
            matrix[k][j] = matrix[k][j] / matrix[k][k];
        }
		for (;j < ROW;j += 4)
		{
			__m128 divee =  _mm_loadu_ps(&matrix[k][j]);
			divee = _mm_div_ps(divee, diver);
			_mm_storeu_ps(&matrix[k][j], divee);
		}
		matrix[k][k] = 1.0;
		for (int i = k + 1; i < ROW; i += 1)
		{
			__m128 mult1 = _mm_load_ps1(&matrix[i][k]);
			int j;
			for (j = k + 1;j < ROW && ((ROW - j) & 3);++j)
            {
                matrix[i][j] = matrix[i][j] - matrix[i][k] * matrix[k][j];
            }
			for (;j < ROW;j += 4)
			{
				__m128 sub1 =  _mm_loadu_ps(&matrix[i][j]);
				__m128 mult2 =  _mm_loadu_ps(&matrix[k][j]);
				mult2 = _mm_mul_ps(mult1, mult2);
				sub1 = _mm_sub_ps(sub1, mult2);
				_mm_storeu_ps(&matrix[i][j], sub1);
			}
			matrix[i][k] = 0.0;
		}
	}
}

void dynamic_omp_cache()
{
    float tmp1,tmp2;
    int i,j,k;
    #pragma omp parallel num_threads(NUM_THREADS), shared(matrix), private(i,j,k,tmp1,tmp2)
    for(i = 0; i < ROW - 1; i++)
    {
        float tmp1 = matrix[i][i];
        #pragma omp single
		for(j = i + 1; j < ROW; j++)
        {
			matrix[i][j] = matrix[i][j] / tmp1;
		}
		matrix[i][i] = 1;
		#pragma omp for schedule(dynamic)
		for(j = i + 1; j < ROW; j++)
		{
		    float tmp2 = matrix[j][i];
			for(k = i + 1; k < ROW; k++)
			{
				matrix[j][k] = matrix[j][k] - matrix[i][k] * tmp2;
			}
			matrix[j][i] = 0;
		}
	}
}

void static_omp_cache()
{
    float tmp1,tmp2;
    int i,j,k;
    #pragma omp parallel num_threads(NUM_THREADS), shared(matrix), private(i,j,k,tmp1,tmp2)
    for(i = 0; i < ROW - 1; i++)
    {
        float tmp1 = matrix[i][i];
        #pragma omp single
		for(j = i + 1; j < ROW; j++)
        {
			matrix[i][j] = matrix[i][j] / tmp1;
		}
		matrix[i][i] = 1;
		#pragma omp for schedule(static)
		for(j = i + 1; j < ROW; j++)
		{
		    float tmp2 = matrix[j][i];
			for(k = i + 1; k < ROW; k++)
			{
				matrix[j][k] = matrix[j][k] - matrix[i][k] * tmp2;
			}
			matrix[j][i] = 0;
		}
	}
}

void guided_omp_cache()
{
    float tmp1,tmp2;
    int i,j,k;
    #pragma omp parallel num_threads(NUM_THREADS), shared(matrix), private(i,j,k,tmp1,tmp2)
    for(i = 0; i < ROW - 1; i++)
    {
        float tmp1 = matrix[i][i];
        #pragma omp single
		for(j = i + 1; j < ROW; j++)
        {
			matrix[i][j] = matrix[i][j] / tmp1;
		}
		matrix[i][i] = 1;
		#pragma omp for schedule(guided)
		for(j = i + 1; j < ROW; j++)
		{
		    float tmp2 = matrix[j][i];
			for(k = i + 1; k < ROW; k++)
			{
				matrix[j][k] = matrix[j][k] - matrix[i][k] * tmp2;
			}
			matrix[j][i] = 0;
		}
	}
}

void dynamic_omp_all()
{
    float tmp1,tmp2;
    int i,j,k;
    #pragma omp parallel num_threads(NUM_THREADS), shared(matrix), private(i,j,k,tmp1,tmp2)
    for(i = 0; i < ROW - 1; i++)
    {
        float tmp1 = matrix[i][i];
        #pragma omp for schedule(dynamic)
		for(j = i + 1; j < ROW; j++)
        {
			matrix[i][j] = matrix[i][j] / tmp1;
		}
		matrix[i][i] = 1;
		#pragma omp for schedule(dynamic)
		for(j = i + 1; j < ROW; j++)
		{
		    float tmp2 = matrix[j][i];
			for(k = i + 1; k < ROW; k++)
			{
				matrix[j][k] = matrix[j][k] - matrix[i][k] * tmp2;
			}
			matrix[j][i] = 0;
		}
	}
}

void dynamic_omp_all_block()
{
    float tmp1,tmp2;
    int i,j,k;
    int block_size = sqrt(ROW * 1.0 / NUM_THREADS);
    #pragma omp parallel num_threads(NUM_THREADS), shared(matrix), private(i,j,k,tmp1,tmp2)
    for(i = 0; i < ROW - 1; i++)
    {
        float tmp1 = matrix[i][i];
        #pragma omp for schedule(dynamic,block_size)
		for(j = i + 1; j < ROW; j++)
        {
			matrix[i][j] = matrix[i][j] / tmp1;
		}
		matrix[i][i] = 1;
		#pragma omp for schedule(dynamic,block_size)
		for(j = i + 1; j < ROW; j++)
		{
		    float tmp2 = matrix[j][i];
			for(k = i + 1; k < ROW; k++)
			{
				matrix[j][k] = matrix[j][k] - matrix[i][k] * tmp2;
			}
			matrix[j][i] = 0;
		}
	}
}

void dynamic_omp_SIMD()
{
    __m128 diver,divee,mult1,mult2,sub1;
    int i,j,k;
    #pragma omp parallel num_threads(NUM_THREADS), shared(matrix), private(i,j,k,diver,divee,mult1,mult2,sub1)
    for(k = 0; k < ROW; ++k)
	{
		diver = _mm_load_ps1(&matrix[k][k]);
		#pragma omp single
		for (j = k + 1;j < ROW && ((ROW - j) & 3);++j)
        {
            matrix[k][j] = matrix[k][j] / matrix[k][k];
        }
		for (;j < ROW;j += 4)
		{
			divee =  _mm_loadu_ps(&matrix[k][j]);
			divee = _mm_div_ps(divee, diver);
			_mm_storeu_ps(&matrix[k][j], divee);
		}
		#pragma omp barrier
		matrix[k][k] = 1.0;
		#pragma omp for schedule(dynamic)
		for (i = k + 1; i < ROW; i += 1)
		{
			mult1 = _mm_load_ps1(&matrix[i][k]);
			int j;
			for (j = k + 1;j < ROW && ((ROW - j) & 3);++j)
            {
                matrix[i][j] = matrix[i][j] - matrix[i][k] * matrix[k][j];
            }
			for (;j < ROW;j += 4)
			{
				sub1 =  _mm_loadu_ps(&matrix[i][j]);
				mult2 =  _mm_loadu_ps(&matrix[k][j]);
				mult2 = _mm_mul_ps(mult1, mult2);
				sub1 = _mm_sub_ps(sub1, mult2);
				_mm_storeu_ps(&matrix[i][j], sub1);
			}
			matrix[i][k] = 0.0;
		}
		#pragma omp barrier
	}
}

void static_omp_SIMD()
{
    __m128 diver,divee,mult1,mult2,sub1;
    int i,j,k;
    #pragma omp parallel num_threads(NUM_THREADS), shared(matrix), private(i,j,k,diver,divee,mult1,mult2,sub1)
    for(k = 0; k < ROW; ++k)
	{
		diver = _mm_load_ps1(&matrix[k][k]);
		#pragma omp single
		for (j = k + 1;j < ROW && ((ROW - j) & 3);++j)
        {
            matrix[k][j] = matrix[k][j] / matrix[k][k];
        }
		for (;j < ROW;j += 4)
		{
			divee =  _mm_loadu_ps(&matrix[k][j]);
			divee = _mm_div_ps(divee, diver);
			_mm_storeu_ps(&matrix[k][j], divee);
		}
		#pragma omp barrier
		matrix[k][k] = 1.0;
		#pragma omp for schedule(static)
		for (i = k + 1; i < ROW; i += 1)
		{
			mult1 = _mm_load_ps1(&matrix[i][k]);
			int j;
			for (j = k + 1;j < ROW && ((ROW - j) & 3);++j)
            {
                matrix[i][j] = matrix[i][j] - matrix[i][k] * matrix[k][j];
            }
			for (;j < ROW;j += 4)
			{
				sub1 =  _mm_loadu_ps(&matrix[i][j]);
				mult2 =  _mm_loadu_ps(&matrix[k][j]);
				mult2 = _mm_mul_ps(mult1, mult2);
				sub1 = _mm_sub_ps(sub1, mult2);
				_mm_storeu_ps(&matrix[i][j], sub1);
			}
			matrix[i][k] = 0.0;
		}
		#pragma omp barrier
	}
}

void guided_omp_SIMD()
{
    __m128 diver,divee,mult1,mult2,sub1;
    int i,j,k;
    #pragma omp parallel num_threads(NUM_THREADS), shared(matrix), private(i,j,k,diver,divee,mult1,mult2,sub1)
    for(k = 0; k < ROW; ++k)
	{
		diver = _mm_load_ps1(&matrix[k][k]);
		#pragma omp single
		for (j = k + 1;j < ROW && ((ROW - j) & 3);++j)
        {
            matrix[k][j] = matrix[k][j] / matrix[k][k];
        }
		for (;j < ROW;j += 4)
		{
			divee =  _mm_loadu_ps(&matrix[k][j]);
			divee = _mm_div_ps(divee, diver);
			_mm_storeu_ps(&matrix[k][j], divee);
		}
		#pragma omp barrier
		matrix[k][k] = 1.0;
		#pragma omp for schedule(guided)
		for (i = k + 1; i < ROW; i += 1)
		{
			mult1 = _mm_load_ps1(&matrix[i][k]);
			int j;
			for (j = k + 1;j < ROW && ((ROW - j) & 3);++j)
            {
                matrix[i][j] = matrix[i][j] - matrix[i][k] * matrix[k][j];
            }
			for (;j < ROW;j += 4)
			{
				sub1 =  _mm_loadu_ps(&matrix[i][j]);
				mult2 =  _mm_loadu_ps(&matrix[k][j]);
				mult2 = _mm_mul_ps(mult1, mult2);
				sub1 = _mm_sub_ps(sub1, mult2);
				_mm_storeu_ps(&matrix[i][j], sub1);
			}
			matrix[i][k] = 0.0;
		}
		#pragma omp barrier
	}
}

void dynamic_omp_SIMD_self()
{
    float tmp1,tmp2;
    int i,j,k;
    #pragma omp parallel num_threads(NUM_THREADS), shared(matrix), private(i,j,k,tmp1,tmp2)
    for(i = 0; i < ROW - 1; i++)
    {
        tmp1 = matrix[i][i];
        #pragma omp single
		for(j = i + 1; j < ROW; j++)
        {
			matrix[i][j] = matrix[i][j] / tmp1;
		}
		matrix[i][i] = 1;
		#pragma omp for schedule(dynamic)
		for(j = i + 1; j < ROW; j++)
		{
		    tmp2 = matrix[j][i];
		    #pragma omp simd
			for(k = i + 1; k < ROW; k++)
			{
				matrix[j][k] = matrix[j][k] - matrix[i][k] * tmp2;
			}
			matrix[j][i] = 0;
		}
	}
}

void gettime(void (*func)())
{
	ll head, tail, freq;
	double time = 0;
	int counter = 0;
	while (INTERVAL > time)
	{
        init();
        QueryPerformanceFrequency((LARGE_INTEGER*)&freq);
        QueryPerformanceCounter((LARGE_INTEGER*)&head);
		func();
		QueryPerformanceCounter((LARGE_INTEGER*)&tail);
		counter++;
		time += (tail - head) * 1000.0 / freq;
	}
	std::cout << time / counter << '\n';
}

int main()
{
	cout<<"plain: ";
	gettime(plain);
	for (NUM_THREADS = 2;NUM_THREADS <= 8;NUM_THREADS++)
	{
		cout<<"threads_num: "<<NUM_THREADS<<endl;
		cout<<"dynamic_omp_cache: ";
		gettime(dynamic_omp_cache);
		cout<<"static_omp_cache: ";
		gettime(static_omp_cache);
		cout<<"guided_omp_cache: ";
		gettime(guided_omp_cache);
		cout<<"dynamic_omp_all: ";
		gettime(dynamic_omp_all);
		cout<<"dynamic_omp_all_block: ";
		gettime(dynamic_omp_all_block);
		cout<<"dynamic_omp_SIMD: ";
		gettime(dynamic_omp_SIMD);
		cout<<"static_omp_SIMD: ";
		gettime(static_omp_SIMD);
		cout<<"guided_omp_SIMD: ";
		gettime(guided_omp_SIMD);
		cout<<"dynamic_omp_SIMD_self:";
		gettime(dynamic_omp_SIMD_self);
	}
	return 0;
}
